Modeling of octave-spanning Kerr frequency combs using a generalized 

mean-field Lugiato-Lefever model 



.1 



Stephane Coen, 1 * Hamish G. Randle, 1 Thibaut Sylvestre, 2 and Miro Erkintalo 1 

1 Physics Department, The University of Auckland, Private Bag 92019, Auckland, New Zealand 
2 Institut FEMTO-ST, Departement d'Optique P. M. Duffieux, Universite de Franche-Comte, CNRS, Besancon, France 

* Corresponding author: s.coen® auckland.ac.nz 



(N 

o 

(N 

> 
O 

ov 

(N 



C/5 

O 

> 

q. 

X. 

Of 



(N- 

>: 

OV 
VO 1 



(N 



X 



Compiled November 30, 2012 

A generalized Lugiato-Lefever equation is numerically solved with a Newton-Raphson method to model Kerr frequency combs. 
We obtain excellent agreement with past experiments, even for an octave-spanning comb. Simulations are much faster than 
with any other technique despite including more modes than ever before. Our study reveals that Kerr combs are associated with 
temporal cavity solitons and dispersive waves, and opens up new avenues for the understanding of Kerr comb formation. © 
2012 Optical Society of America 

OCIS codes: 230.5750, 190.5530, 190.4380 



First observed in 2007, frequency comb generation in mono- 
lithic microring resonators has aroused significant research 
interest Q]|2|- A minuscule footprint, power efficiency, and 
CMOS -compatibility make said structures ideal for on-chip 
frequency comb generation. Applications range from spec- 
troscopy to telecommunications. 

Although comb generation in high-Q Kerr resonators has 
been extensively reported experimentally, theoretical studies 
are comparatively scarce. In part this deficiency can be linked 
to the intractable computational complexity of the existing 
models. Indeed, the use of a nonlinear Schrodinger (NLS) 
equation and appropriate coupling-mediated boundary condi- 
tions requires hundreds of millions of roundtrips for conver- 
gence, owing to the exceedingly high Q of the structures J3|. 
Likewise the number of four- wave mixing-mediated coupling 
terms in coupled mode equation models scales cubically with 
the number of modes, limiting such modeling to narrowband 
combs [4|. Matsko et al. Q also considered an analytic ap- 
proach but it is restricted to resonators of infinite intrinsic 
Q-factor and with group-velocity dispersion (GVD) limited 
to 2nd order. Clearly, a growing demand exists for realistic 
yet computationally efficient methods for the modeling of fre- 
quency combs in high-Q resonators. 

In this Letter we report on direct numerical modeling of 
Kerr frequency combs using a generalized Lugiato-Lefever 
(LL) equation J6), and find good agreement with reported ex- 
periments. Significantly, the conducted computations are far 
less intensive than previous methods, allowing for the rapid 
modeling of octave-spanning combs with arbitrarily low rep- 
etition rates on a consumer grade computer. We also show 
that the results obtained from the proposed model provide 
significant insights into Kerr combs. In particular, we high- 
light how localized dissipative structures known as temporal 
cavity solitons (CSs) [7] are stable stationary solutions of the 
governing equation, and how specific comb features can be 
intuitively described in terms of CS propagation. We expect 
the reported technique to become an invaluable tool for the 
understanding and tailoring of nonlinear optical phenomena 
in high-Q resonators. 



We consider a typical ring resonator configuration (Fig.Q~|i: 
a continuous-wave (cw) driving field E- m with power P m = 
|£in| 2 is coherently added to the lightwave circulating in the 
resonator through a coupler with power transmission coeffi- 
cient 9. The fourth port of the coupler is used to extract the 
output field. Mathematically the intracavity field £'( m+1 ) (0, t) 
at the beginning of the (m + 1 ) th roundtrip can be related to 
the field (L, t) at the end of the m th roundtrip as: 



£ (m+1) (0,T) = V9E in + VT~eE^ m \L, T)< 



Mi 



(1) 



where i is the time, L the roundtrip-length of the resonator, 
and 0o the linear phase accumulated by the intracavity field 
with respect to the the pump field over one roundtrip. 

Assuming that light propagates in a single spatial mode, the 
evolution of the field through one roundtrip of the resonator 
is governed by the well-known (generalized) NLS equation: 



dE(z,r) _ Oj 



dz 



2 tk k\ \ dz 
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E + iy\E\ 2 E. (2) 



Here (X{ is the linear absorption coefficient inside the res- 
onator, ^ = d k p /d(O k \ a=(iK) the dispersion coefficients as- 
sociated with the Taylor series expansion of the propagation 
constant J3 (co) about the center frequency C0q of the driving 
field, and y = 7j?COo/(cA e ff) the nonlinearity coefficient with 
«2 the nonlinear refractive index and A e g the effective modal 
area of the resonator mode. 

The boundary conditions Eq. (Q]) combined with the NLS 
Eq. (O form an infinite-dimensional map that describes com- 
pletely the dynamics of a ring resonator of any size or shape 
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Fig. 1. (Color online) Schematic of the resonator configuration. 
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(toroid, racetrack, . . .)■ In low loss structures, the field varies 
only slightly between consecutive roundtrips, making direct 
simulations of these equations very slow yj. In these condi- 
tions, however, it is well known that Eqs. (H)-© can be aver- 
aged into an externally driven NLS equation (see, e.g., J8)), 



(a) Simulation 



(b) Experiment 
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where t^ is the roundtrip time, a = (cCi + 9)/2 describes the 
total cavity losses, and 5q = 2nl — 0o is the detuning with / the 
order of the cavity resonance closest to the driving field. The 
continuous variable t measures the slow time of the cavity, 
and can be linked to the roundtrip index as E(t = mtR,T) — 
£( m )(0,T). Equation (0 coincides with the master equation 
of and, with fa = for k > 3, is formally identical to the 
mean-field LL model of a diffractive cavity l6ll9l [T0i . It has 
also been extensively used for the description of passive cav- 
ities constructed of sing le-mode fibers f7ll8 TTT1fl2l . In partic- 
ular, spatial pattern formation and the so-called modulation 
instability (MI) studied in some of these earlier works can be 
directly connected to frequency comb generation 16191121 . MI 
was also shown to occur in the normal GVD regime in pres- 
ence of cavity boundary conditions I0QT), wn i cn is directly 
relevant to combs. We must also note that the expressions of 
the characteristic bistable response of the LL model as well 
as that of the intracavity MI gain [8) in fact coincide after 
normalization with corresponding results obtained from the 
coupled-mode equations of JH, suggesting an intrinsic link 
between the two approaches. 

The generalized LL Eq. (0 holds in the limit of high fi- 
nesse cavities & 1. For typical high-Q resonators, the fi- 
nesse & ~ 10 2 -10 5 . Also, dispersion must be "weak" over 
one roundtrip, Y,k>2 faLA(O k /k\ < n, where A© is the (an- 
gular) spectral width of the generated comb. This condition 
was found to be satisfied a posteriori for all the results dis- 
cussed below, thereby asserting the validity of Eq. (0 for the 
description of Kerr combs. The LL equation has two substan- 
tial advantages compared to the infinite dimensional map Eqs. 
(Q]l-©- On the one hand, Eq. (0 can be numerically inte- 
grated with the split-step Fourier method using an integra- 
tion step corresponding to multiple roundtrips, significantly 
reducing the computational burden in obtaining steady state 
solutions. On the other hand, the steady state solutions can be 
obtained directly by setting dE /dt = and looking for the 
roots of the right-hand-side of Eq. (0. Although the latter ap- 
proach does not yield insights into the dynamics of comb for- 
mation, it is computationally orders of magnitude more effi- 
cient than split-step integration. Here we restrict our attention 
to stationary solutions obtained by a multidimensional root- 
finding Newton-Rhapson method. Derivatives are computed 
with Fourier transforms and the span of the temporal window 
coincides with meaning the samples of our frequency grid 
are spaced by the free-spectral range, FSR = 1 /?r. 

As a first example, we plot in Fig. 0a) the intracavity spec- 
trum of a stable steady-state solution of Eq. (0 obtained with 
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Fig. 2. (a) Steady state solution of Eq. lO for a critically-coupled, 
3.8 mm diameter MgF2 whispering gallery mode resonator with a 
40 [im mode-field diameter and a loaded Q = 1.90 ■ 10 9 . FSR = 



18.2 GHz; y = 0.032 W^km -1 ; j3 2 



13 ps km 
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1.75- 10~ 5 ; P m = 55.6 mW; L = 11.9 mm; Sq = 0.0012. (b) Corre- 
sponding experimental spectrum after jl3| . 

simulation parameters listed in the caption and approximating 
the experimental values of ifTJl . We note that some of these 
parameters have large uncertainties but this does not invali- 
date our conclusions. Figure0b) is the corresponding exper- 
imental output spectrum, and clearly excellent agreement is 
observed (abstraction must be made of the pump mode inten- 
sity which is modified at the output by the reflected pump). 
Some of the discrepancies could be traced to effects unac- 
counted for in Eq. (0 such as wavelength-dependent losses 
or overlap integrals, or to experimental fluctuations. We also 
show the temporal profile of the solution as an inset of Fig. 
a) so as to highlight that the intracavity field corresponds to 
a ~ 400 fs-duration pulse with a peak power of 100 mW atop 
a weak cw background. It is well-known that the LL equation 
possesses such localized-CS solutions repeating at the cavity 
repetition rate thus forming a frequency comb in the spectral 
domain J7][10). In fact, in solving Eq. (0, we have not found 
any type of stable steady-state comb solutions not made up 
of single or multiple CSs. This strongly suggests that stable 
frequency combs generated in high-Q cavities generally cor- 
respond to trains of CSs which is in good accordance with 
recent studies [00Q3) and was already suggested in |7). Our 
analysis also often reveals coexisting unstable states associ- 
ated with MI and breathing CSs which can in practice pre- 
clude the observation of stable combs. Care must therefore 
be taken in interpreting some reported experimental results: 
They may be associated with rapid fluctuations and only ap- 
pear stable because of the averaging of spectrum analyzers. 

Our next example (Fig. is an octave-spanning comb, 
strongly influenced by higher-order dispersion. The simula- 
tion parameters are radically different and are approximated 
from the experiment of IPT31 (see caption). Again, based on 
the same model Eq. (0, we observe an excellent agreement 
between the stable steady-state numerical spectrum plotted 
in Fig. 0a) and the experimental measurements (see Fig. 2 
in 03]). At this point we must emphasize that, despite the 
large number of spectral modes included (1024), the compu- 
tation time in obtaining the results shown both in Figs. 0a) 
and 0a) was of the order of seconds on a standard computer. 
It is precisely this unprecedented speed that is the key ad- 
vantage of our technique. In fact, to our knowledge, the Kerr 
comb of Fig. 0a) has the largest number of optical modes ob- 
tained from a theoretical model. Simulating octave-spanning 
combs with sub-40 GHz repetiton rates would necessitate to 
ramp up the number of modes to 4096. In this case, computa- 
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Fig. 3. (Color online) (a) Steady-state solution of Eq. (|3j with 
parameters mimicking a critically-coupled, 200 flm diameter sili- 
con nitride resonator with a loaded Q = 3- 10 5 approximated from 
fBl . Dispersion as per (b); FSR = 226 GHz; y = 1 W^'m" 1 ; 
a = 6 = 0.009; P m = 755 mW; L = 628 jlim; 5q = 0.0534. (c) Time- 
frequency representation of the simulation result calculated using a 
100 fs gate function. Subsequent roundtrips are separated by verti- 
cal lines whilst the horizontal line indicates the predicted Cerenkov 
wavelength. 

tion time increases to about 2-3 minutes (or less if consider- 
ing neighboring solutions), but is still very reasonable. 

A particular feature seen in Fig.[3ja) is the strong narrow- 
band low-frequency component centered about 2150 nm. A 
similar feature can be witnessed in the corresponding experi- 
mental measurements ifTBI . and also in other previous experi- 
ments (see, e.g., Fig. 1 in [16]). We interpret this component 
as a Cerenkov-like resonant dispersive wave (DW) emitted 
by the CS circulating in the resonator ifTTl . To show this ex- 
plicitly, we plot in Fig. 0c) the time-frequency representation 
of the intracavity electric field. Here we exploit the periodic 
boundary conditions, and expand the fast temporal axis across 
three cavity roundtrips. We can clearly see how the intracavity 
field consists of a train of pulses atop a cw background (i.e., 
CSs) and we can identify the narrowband 2150 nm compo- 
nent as DWs emitted by individual CSs. This is further high- 
lighted by the dashed horizontal line in Fig. 0c) which in- 
dicates the predicted wavelength Adw of a Cerenkov-wave. 
Specifically, neglecting the nonlinear contribution, the perti- 
nent phase-matching condition governing the resonant DW 
is J3(c0d W ) = /3(cocs) -j3i(cocs) • (e>cs- g>dw) 03, where 
©dw = 27Tc/A,dw and <x>cs are the central (angular) frequen- 
cies of the DW and the CS, respectively. With our numerical 
parameters, this condition yields Adw = 2149 nm, in excel- 
lent agreement with the observed spectral peak in Fig. 0a). 
Because a (cavity) soliton in the anomalous GVD regime can 
only excite DWs in the normal GVD regime, our observations 
suggest that the ubiquitous asymmetry of Kerr combs towards 
the normal GVD regime may in fact be explained by the exci- 
tation of resonant DWs by CSs in the anomalous GVD region 
in a way akin to supercontinuum generation 11711 . Finally, we 
note that DW emission has recently been described in terms 
of cascaded four-wave mixing, which further highlights its 
relevance to Kerr combs lfl8l . 

In conclusion, we have used a generalized LL equation 



to model high-Q resonator Kerr frequency combs. It can be 
solved with a Newton-Raphson solver, providing results in 
a matter of seconds, much faster than any other technique, 
and for widely different cases, all the way to octave-spanning 
combs with hundreds of modes. Our results are in good agree- 
ment with experiments, proving that the LL model captures 
the essential physics of Kerr combs. We also find that Kerr 
frequency combs are associated in the temporal domain with 
CSs J7| and associated DWs. CSs are well known localized 
dissipative structures which have been studied extensively in 
the past, in particular in the spatial domain Il9l [l0ll . so we 
anticipate that the LL model will provide both a deeper un- 
derstanding to Kerr frequency combs and computational effi- 
ciency. The Newton method also provides information about 
dynamical instabilities of Kerr combs through an eigenvalue 
analysis of the Jacobian of the system, which is obtained at 
no extra cost. Some of our preliminary results highlight the 
presence of Hopf bifurcations associated with breathing CSs, 
i.e., combs with periodically modulated features. 

We thank Dr I. S. Grudinin for kindly providing experi- 
mental data. We also acknowledge support from the Marsden 
Fund of The Royal Society of New Zeland. 
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